home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / prfalign.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  24.6 KB  |  949 lines  |  [TEXT/R*ch]
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <math.h>
  5. #include "clustalw.h"
  6. #define ENDALN 127
  7.  
  8. #define MAX(a,b) ((a)>(b)?(a):(b))
  9. #define MIN(a,b) ((a)<(b)?(a):(b))
  10.  
  11. /*
  12.  *   Prototypes
  13.  */
  14.  
  15. extern void calc_prf1(int **Profile, char **alignment, int *gaps, 
  16.             int matrix[NUMRES][NUMRES],
  17.              int *weight, int prf_length, int ThisSeq, int LastSeq);
  18. extern void calc_prf2(int **Profile, char **alignment, int *weight,
  19.                       int prf_length, int ThisSeq, int LastSeq);
  20. extern void calc_gap_coeff(char **alignment, int *gaps, int **Profile,
  21.                            int FirstSeq, int LastSeq, int prf_length,
  22.                            int gapcoef, int lencoef);
  23. extern int  get_matrix(int *matptr, int *xref, int matrix[NUMRES][NUMRES], int neg_flag);
  24.  
  25.  
  26. extern void *ckalloc(size_t bytes);
  27. extern void *ckrealloc(void *ptr, size_t bytes);
  28. extern void ckfree(void *);
  29.  
  30. int prfalign(int *group, int *aligned);
  31. static int     pdiff(short A,short B,short i,short j);
  32. static int     prfscore(short n, short m);
  33. static int     gap_penalty1(short i, short j,short k);
  34. static int     open_penalty1(short i, short j);
  35. static int     ext_penalty1(short i, short j);
  36. static int     gap_penalty2(short i, short j,short k);
  37. static int     open_penalty2(short i, short j);
  38. static int     ext_penalty2(short i, short j);
  39. static void     padd(short k);
  40. static void     pdel(short k);
  41. static void     palign(void);
  42. static void     ptracepath(short se1, short se2, short *alen);
  43. static void     add_ggaps(int n1, int n2);
  44.  
  45. /*
  46.  *   Global variables
  47.  */
  48. extern double     **tmat;
  49. extern float     gap_open, gap_extend;
  50. extern int     gap_pos1, gap_pos2, max_aa;
  51. extern int     nseqs;
  52. extern int     *seqlen_array;
  53. extern int     *seq_weight;
  54. extern int        debug;
  55. extern int     neg_matrix;
  56. extern int     mat_avscore;
  57. extern int      blosum30mt[], blosum35mt[], blosum40mt[], blosum45mt[];
  58. extern int      blosum50mt[], blosum55mt[], blosum62mt2[], blosum70mt[];
  59. extern int      blosum75mt[], blosum80mt[], blosum85mt[], blosum90mt[];
  60. extern int      pam20mt[], pam60mt[];
  61. extern int      pam120mt[], pam160mt[], pam250mt[], pam350mt[];
  62. extern int      md_350mt[], md_250mt[], md_120mt[], md_40mt[];
  63. extern int      idmat[], usermat[];
  64. extern int    def_aa_xref[], aa_xref[];
  65. extern Boolean    distance_tree;
  66. extern Boolean    dnaflag, is_weight;
  67. extern char     mtrxname[];
  68. extern char     **seq_array;
  69. extern char     *amino_acid_codes;
  70.  
  71. static int     print_ptr,last_print;
  72. /* static int     displ[2*MAXLEN+1]; */
  73. extern int   *displ;
  74.  
  75.  
  76. static char       **alignment;
  77. static int        *aln_len;
  78. static int        *aln_weight;
  79. static char       *aln_path1, *aln_path2;
  80. static short        alignment_len;
  81. static int        **profile1, **profile2;
  82. static int        matrix[NUMRES][NUMRES];
  83. static int        nseqs1, nseqs2;
  84. static int        prf_length1, prf_length2;
  85. static int        verbose = FALSE;
  86. /* static int        gaps[MAXLEN]; */
  87. extern int * gaps;
  88.  
  89. static int        gapcoef;
  90. static int        lencoef;
  91.  
  92. int prfalign(int *group, int *aligned)
  93. {
  94.  
  95.   static int    i, j, k, count = 0;
  96.   static int    NumSeq,temp,len,insert;
  97.   static int    winlength;
  98.   static short    se1, se2, sb1, sb2;
  99.   static int    maxres;
  100.   static int    *matptr;
  101.   static int    *mat_xref;
  102.   static char   c;
  103.   static char   reply[FILENAMELEN+1];
  104.   static int    score;
  105.   static float  scale,ftemp;
  106.   static double logmin,logdiff;
  107.   static double pcid;
  108.  
  109.  
  110.   alignment = (char **) ckalloc( nseqs * sizeof (char *) );
  111.   aln_len = (int *) ckalloc( nseqs * sizeof (int) );
  112.   aln_weight = (int *) ckalloc( nseqs * sizeof (int) );
  113.  
  114.   for (i=0;i<nseqs;i++)
  115.      if (aligned[i+1] == 0) group[i+1] = 0;
  116.  
  117.   nseqs1 = nseqs2 = 0;
  118.   for (i=0;i<nseqs;i++)
  119.     {
  120.         if (group[i+1] == 1) nseqs1++;
  121.         else if (group[i+1] == 2) nseqs2++;
  122.     }
  123.  
  124.   if ((nseqs1 == 0) || (nseqs2 == 0)) return(0.0);
  125.  
  126.   if (nseqs2 > nseqs1)
  127.     {
  128.      for (i=0;i<nseqs;i++)
  129.        {
  130.           if (group[i+1] == 1) group[i+1] = 2;
  131.           else if (group[i+1] == 2) group[i+1] = 1;
  132.        }
  133.     }
  134.  
  135.   if (dnaflag)
  136.      {
  137.        for(i=0;i<5;++i)
  138.          {
  139.             for(j=0;j<5;++j) {
  140.                  if(i==j)
  141.                      matrix[i][j]=1000;
  142.                  else
  143.                      matrix[j][i]=0;
  144.              }
  145.             if(is_weight) {
  146.                 matrix[1][3]=500;
  147.                 matrix[1][4]=500;
  148.                 matrix[3][1]=500;
  149.                 matrix[4][1]=500;
  150.                 matrix[2][0]=500;
  151.                 matrix[0][2]=500;
  152.              }
  153.          }
  154.     }
  155.   else
  156.     {
  157.  
  158. /*
  159.    calculate the mean of the sequence pc identities between the two groups
  160. */
  161.         count = 0;
  162.         pcid = 0.0;
  163.         for (i=0;i<nseqs;i++)
  164.           {
  165.              if (group[i+1] == 1)
  166.              for (j=0;j<nseqs;j++)
  167.                if (group[j+1] == 2)
  168.                     {
  169.                        count++;
  170.                        if (pcid < tmat[i+1][j+1]) pcid = tmat[i+1][j+1];
  171. /*
  172.                        pcid += tmat[i+1][j+1];
  173. */
  174.                     }
  175.           }
  176. /*
  177.   pcid = pcid/(float)count;
  178. */
  179.  
  180. if (debug>0) fprintf(stderr,"mean tmat %3.1f\n", pcid);
  181.  
  182.        scale = 0.5;
  183.        if (strcmp(mtrxname, "blosum") == 0)
  184.         {
  185.            if (distance_tree == FALSE) matptr = blosum40mt;
  186.            else if (pcid > 80.0)
  187.              {
  188.                 matptr = blosum80mt;
  189.                 scale *= -3.5 + pcid/20.0;
  190.              }
  191.            else if (pcid > 60.0)
  192.              {
  193.                 matptr = blosum62mt2;
  194.                 scale *= -2.5 + pcid/20.0;
  195.              }
  196.            else if (pcid > 30.0)
  197.              {
  198.                 matptr = blosum45mt;
  199.                 scale *= -0.4 + pcid/30.0;
  200.              }
  201.            else if (pcid > 10.0)
  202.              {
  203.                 matptr = blosum30mt;
  204.                 scale *= pcid/30.0;
  205.              }
  206.            else
  207.              {
  208.                 matptr = blosum30mt;
  209.                 scale *= 0.30;
  210.              }
  211.            mat_xref = def_aa_xref;
  212.  
  213.         }
  214.        else if (strcmp(mtrxname, "md") == 0)
  215.         {
  216.            if (distance_tree == FALSE) matptr = md_250mt;
  217.            else if (pcid > 80.0) matptr = md_40mt;
  218.            else if (pcid > 60.0) matptr = md_120mt;
  219.            else if (pcid > 40.0) matptr = md_250mt;
  220.            else matptr = md_350mt;
  221.            mat_xref = def_aa_xref;
  222.         }
  223.        else if (strcmp(mtrxname, "pam") == 0)
  224.         {
  225.            if (distance_tree == FALSE) matptr = pam120mt;
  226.            else if (pcid > 80.0) matptr = pam20mt;
  227.            else if (pcid > 60.0) matptr = pam60mt;
  228.            else if (pcid > 40.0) matptr = pam120mt;
  229.            else matptr = pam350mt;
  230.            mat_xref = def_aa_xref;
  231.         }
  232.        else if (strcmp(mtrxname, "id") == 0)
  233.         {
  234.            matptr = idmat;
  235.            mat_xref = def_aa_xref;
  236.         }
  237.        else 
  238.         {
  239.            matptr = usermat;
  240.            mat_xref = aa_xref;
  241.         }
  242.  
  243.       maxres = get_matrix(matptr, mat_xref, matrix, neg_matrix);
  244.       if (maxres == 0)
  245.         {
  246.            fprintf(stderr,"Error: matrix %s not found\n", mtrxname);
  247.            return(-1.0);
  248.         }
  249.     }
  250.  
  251. /*
  252.   Make the first profile.
  253. */
  254.   prf_length1 = 0;
  255.   nseqs1 = 0;
  256.   for (i=0;i<nseqs;i++)
  257.     {
  258.        if (group[i+1] == 1)
  259.           {
  260.              len = seqlen_array[i+1];
  261.              alignment[nseqs1] = (char *) ckalloc( (len+2) * sizeof (char) );
  262.              for (j=0;j<len;j++)
  263.                alignment[nseqs1][j] = seq_array[i+1][j+1];
  264.              alignment[nseqs1][j] = ENDALN;
  265.              aln_len[nseqs1] = len;
  266.              aln_weight[nseqs1] = seq_weight[i];
  267.              if (len > prf_length1) prf_length1 = len;
  268.              nseqs1++;
  269.           }
  270.     }
  271.  
  272. /*
  273.   Make the second profile.
  274. */
  275.   prf_length2 = 0;
  276.   nseqs2 = 0;
  277.   for (i=0;i<nseqs;i++)
  278.     {
  279.        if (group[i+1] == 2)
  280.           {
  281.              len = seqlen_array[i+1];
  282.              alignment[nseqs1+nseqs2] =
  283.                    (char *) ckalloc( (len+2) * sizeof (char) );
  284.              for (j=0;j<len;j++)
  285.                alignment[nseqs1+nseqs2][j] = seq_array[i+1][j+1];
  286.              alignment[nseqs1+nseqs2][j] = ENDALN;
  287.              aln_len[nseqs1+nseqs2] = len;
  288.              aln_weight[nseqs1+nseqs2] = seq_weight[i];
  289.              if (len > prf_length2) prf_length2 = len;
  290.              nseqs2++;
  291.           }
  292.     }
  293.  
  294. if (debug>0) fprintf(stderr,"average: %d\n",(pint)mat_avscore);
  295.  
  296.   logmin = log((double)(MIN(prf_length1,prf_length2)));
  297.   if (prf_length1<=prf_length2)
  298.      logdiff = 1.0-log((double)((float)prf_length1/(float)prf_length2));
  299.   else
  300.      logdiff = 1.0-log((double)((float)prf_length2/(float)prf_length1));
  301.  
  302. /*
  303.     round logdiff to the nearest integer
  304. */
  305.   if ((logdiff-(int)logdiff) > 0.5) logdiff=ceil(logdiff);
  306.   else                             logdiff=floor(logdiff);
  307.  
  308. if (debug>0) fprintf(stderr,"%d %d logmin %f   logdiff %f\n",
  309. (pint)prf_length1,(pint)prf_length2, logmin,logdiff);
  310.  
  311.   if (dnaflag)
  312.     {
  313.           gapcoef = 100.0 * gap_open;
  314.           lencoef = 100.0 * gap_extend;
  315.     }
  316.   else
  317.     {
  318.      if (neg_matrix == TRUE)
  319.        {
  320.           if (mat_avscore <= 0)
  321.               gapcoef = 200.0 * (gap_open + logmin);
  322.           else
  323.               gapcoef = 200.0 * (gap_open + logmin);
  324.           lencoef = 200.0 * gap_extend * logdiff;
  325.        }
  326.      else
  327.        {
  328.           if (mat_avscore <= 0)
  329.               gapcoef = 100.0 * (float)(gap_open + logmin);
  330.           else
  331.               gapcoef = scale * mat_avscore * (float)(gap_open + logmin);
  332.           lencoef = 100.0 * gap_extend * logdiff;
  333.        }
  334.     }
  335.  
  336. if (debug>0)
  337. {
  338. fprintf(stderr,"Gap Open %d    Gap Extend %d\n", (pint)gapcoef,(pint)lencoef);
  339. fprintf(stderr,"Window length %d\n", (pint)winlength);
  340. fprintf(stderr,"Matrix  %s\n", mtrxname);
  341. }
  342.  
  343.   profile1 = (int **) ckalloc( (prf_length1+2) * sizeof (int *) );
  344.   for(i=0; i<prf_length1+2; i++)
  345.        profile1[i] = (int *) ckalloc( (LENCOL+2) * sizeof(int) );
  346.  
  347.   profile2 = (int **) ckalloc( (prf_length2+2) * sizeof (int *) );
  348.   for(i=0; i<prf_length2+2; i++)
  349.        profile2[i] = (int *) ckalloc( (LENCOL+2) * sizeof(int) );
  350.  
  351. /*
  352.   calculate the Gap Coefficients.
  353. */
  354.  
  355.      calc_gap_coeff(alignment, gaps, profile1,
  356.            0, nseqs1, prf_length1, gapcoef, lencoef);
  357.  
  358. /*
  359.   calculate the profile matrix.
  360. */
  361.      calc_prf1(profile1, alignment, gaps, matrix,
  362.           aln_weight, prf_length1, 0, nseqs1);
  363.  
  364. if (debug>4)
  365. {
  366. extern char *amino_acid_codes;
  367.   for (j=0;j<=max_aa;j++)
  368.     fprintf(stderr,"%c    ", amino_acid_codes[j]);
  369.  fprintf(stderr,"\n");
  370.   for (i=0;i<prf_length1;i++)
  371.    {
  372.     for (j=0;j<=max_aa;j++)
  373.       fprintf(stderr,"%d ", (pint)profile1[i+1][j]);
  374.     fprintf(stderr,"%d ", (pint)profile1[i+1][gap_pos1]);
  375.     fprintf(stderr,"%d ", (pint)profile1[i+1][gap_pos2]);
  376.     fprintf(stderr,"%d %d\n",(pint)profile1[i+1][GAPCOL],(pint)profile1[i+1][LENCOL]);
  377.    }
  378. }
  379.  
  380. /*
  381.   calculate the Gap Coefficients.
  382. */
  383.  
  384.      calc_gap_coeff(alignment, gaps, profile2,
  385.            nseqs1, nseqs1+nseqs2, prf_length2, gapcoef, lencoef);
  386.  
  387. /*
  388.   calculate the profile matrix.
  389. */
  390.      calc_prf2(profile2, alignment, aln_weight,
  391.            prf_length2, nseqs1, nseqs1+nseqs2);
  392.  
  393.   aln_path1 = (char *) ckalloc( (prf_length1 + prf_length2) * sizeof(char) );
  394.   aln_path2 = (char *) ckalloc( (prf_length1 + prf_length2) * sizeof(char) );
  395.  
  396. if (debug>4)
  397. {
  398. extern char *amino_acid_codes;
  399.   for (j=0;j<=max_aa;j++)
  400.     fprintf(stderr,"%c    ", amino_acid_codes[j]);
  401.  fprintf(stderr,"\n");
  402.   for (i=0;i<prf_length2;i++)
  403.    {
  404.     for (j=0;j<=max_aa;j++)
  405.       fprintf(stderr,"%d ", (pint)profile2[i+1][j]);
  406.     fprintf(stderr,"%d ", (pint)profile2[i+1][gap_pos1]);
  407.     fprintf(stderr,"%d ", (pint)profile2[i+1][gap_pos2]);
  408.     fprintf(stderr,"%d %d\n",(pint)profile2[i+1][GAPCOL],(pint)profile2[i+1][LENCOL]);
  409.    }
  410. }
  411.  
  412. /*
  413.    align the profiles
  414. */
  415. /* use Myers and Miller to align two sequences */
  416.  
  417.   last_print = 0;
  418.   print_ptr = 1;
  419.  
  420.   sb1 = sb2 = 0;
  421.   se1 = prf_length1;
  422.   se2 = prf_length2;
  423.  
  424.   score = pdiff(sb1, sb2, se1-sb1, se2-sb2);
  425.  
  426.   ptracepath(se1, se2, &alignment_len);
  427.  
  428.   add_ggaps(nseqs1, nseqs2);
  429.  
  430.   for (i=0;i<prf_length1+2;i++)
  431.      ckfree((void *)profile1[i]);
  432.   ckfree((void *)profile1);
  433.  
  434.   for (i=0;i<prf_length2+2;i++)
  435.      ckfree((void *)profile2[i]);
  436.   ckfree((void *)profile2);
  437.  
  438.   prf_length1 = alignment_len;
  439.  
  440.   ckfree((void *)aln_path1);
  441.   ckfree((void *)aln_path2);
  442.  
  443.   NumSeq = 0;
  444.   for (j=0;j<nseqs;j++)
  445.     {
  446.        if (group[j+1]  == 1)
  447.          {
  448.             seqlen_array[j+1] = prf_length1;
  449.             for (i=0;i<prf_length1;i++)
  450.               seq_array[j+1][i+1] = alignment[NumSeq][i];
  451.             NumSeq++;
  452.          }
  453.     }
  454.   for (j=0;j<nseqs;j++)
  455.     {
  456.        if (group[j+1]  == 2)
  457.          {
  458.             seqlen_array[j+1] = prf_length1;
  459.             for (i=0;i<prf_length1;i++)
  460.               seq_array[j+1][i+1] = alignment[NumSeq][i];
  461.             NumSeq++;
  462.          }
  463.     }
  464.  
  465.   for (i=0;i<nseqs1+nseqs2;i++)
  466.      ckfree((void *)alignment[i]);
  467.   ckfree((void *)alignment);
  468.  
  469.   ckfree((void *)aln_len);
  470.  
  471.   return(score/100);
  472.  
  473. }
  474.  
  475. static void add_ggaps(int n1, int n2)
  476. {
  477.    int i,j,k,ix;
  478.    int len;
  479.    char *ta;
  480.  
  481.    ta = (char *) ckalloc( MAXLEN * sizeof (char) );
  482.  
  483.    for (j=0;j<nseqs1;j++)
  484.      {
  485.       ix = 0;
  486.       for (i=0;i<alignment_len;i++)
  487.         {
  488.            if (aln_path1[i] == 2)
  489.               {
  490.                  if (ix < aln_len[j])
  491.                     ta[i] = alignment[j][ix];
  492.                  else 
  493.                     ta[i] = ENDALN;
  494.                  ix++;
  495.               }
  496.            else if (aln_path1[i] == 1)
  497.               {
  498. /*
  499.    insertion in first alignment...
  500. */
  501.                  ta[i] = gap_pos1;
  502.               }
  503.            else
  504.               {
  505.                  fprintf(stderr,"Error in aln_path\n");
  506.               }
  507.          }
  508.        ta[i] = ENDALN;
  509.        
  510.        len = alignment_len;
  511.        if (len >= MAXLEN)
  512.       {
  513.              fprintf(stderr,"Error: alignment is too long (Max. %d)\n",MAXLEN);
  514.              exit(1);
  515.           }
  516.        alignment[j] = (char *)ckrealloc(alignment[j], (len+2) * sizeof (char));
  517.        for (i=0;i<len;i++)
  518.          alignment[j][i] = ta[i];
  519.        alignment[j][i] = ENDALN;
  520.        aln_len[j] = len;
  521.       }
  522.  
  523.  
  524.    for (j=nseqs1;j<nseqs1+nseqs2;j++)
  525.      {
  526.       ix = 0;
  527.       for (i=0;i<alignment_len;i++)
  528.         {
  529.            if (aln_path2[i] == 2)
  530.               {
  531.                  if (ix < aln_len[j])
  532.                     ta[i] = alignment[j][ix];
  533.                  else 
  534.                     ta[i] = ENDALN;
  535.                  ix++;
  536.               }
  537.            else if (aln_path2[i] == 1)
  538.               {
  539. /*
  540.    insertion in second alignment...
  541. */
  542.                  ta[i] = gap_pos1;
  543.               }
  544.            else
  545.               {
  546.                  fprintf(stderr,"Error in aln_path\n");
  547.               }
  548.          }
  549.        ta[i] = ENDALN;
  550.        
  551.        ckfree((void *)alignment[j]);
  552.        len = alignment_len;
  553.        alignment[j] = (char *) ckalloc( (len+2) * sizeof (char) );
  554.        for (i=0;i<len;i++)
  555.          alignment[j][i] = ta[i];
  556.        alignment[j][i] = ENDALN;
  557.        aln_len[j] = len;
  558.       }
  559.  
  560.    ckfree((void *)ta);
  561.  
  562. if (debug>0)
  563. {
  564.   char c;
  565.   extern char *amino_acid_codes;
  566.  
  567.    for (i=0;i<nseqs1+nseqs2;i++)
  568.      {
  569.       for (j=0;j<alignment_len;j++)
  570.        {
  571.         if (alignment[i][j] == ENDALN) break;
  572.         else if ((alignment[i][j] < 0) || (alignment[i][j] > max_aa))  c = '-';
  573.         else c = amino_acid_codes[alignment[i][j]];
  574.         fprintf(stderr,"%c", c);
  575.        }
  576.       fprintf(stderr,"\n\n");
  577.      }
  578. }
  579.  
  580. }                  
  581.  
  582. static int prfscore(short n, short m)
  583. {
  584.    short    ix, iy;
  585.    int  score;
  586.  
  587.    score = 0.0;
  588.    for (ix=0; ix<=max_aa; ix++)
  589.      {
  590.          score += (profile1[n][ix] * profile2[m][ix])/10;
  591.      }
  592.    score += (profile1[n][gap_pos1] * profile2[m][gap_pos1])/10;
  593.    score += (profile1[n][gap_pos2] * profile2[m][gap_pos2])/10;
  594.    return(score);
  595.    
  596. }
  597.  
  598. static void ptracepath(short se1, short se2, short *alen)
  599. {
  600.   short i,j,k,i1,i2,pos,to_do,len;
  601.  
  602.     pos = 0;
  603.  
  604.     to_do=print_ptr-1;
  605.  
  606.     for(i=1;i<=to_do;++i) {
  607. if (debug>0) fprintf(stderr,"%d ",(pint)displ[i]);
  608.             if(displ[i]==0) {
  609.                     aln_path1[pos]=2;
  610.                     aln_path2[pos]=2;
  611.                     ++pos;
  612.             }
  613.             else {
  614.                     if((k=displ[i])>0) {
  615.                             for(j=0;j<=k-1;++j) {
  616.                                     aln_path2[pos+j]=2;
  617.                                     aln_path1[pos+j]=1;
  618.                             }
  619.                             pos += k;
  620.                     }
  621.                     else {
  622.                             k = (displ[i]<0) ? displ[i] * -1 : displ[i];
  623.                             for(j=0;j<=k-1;++j) {
  624.                                     aln_path1[pos+j]=2;
  625.                                     aln_path2[pos+j]=1;
  626.                             }
  627.                             pos += k;
  628.                     }
  629.             }
  630.     }
  631. if (debug>0) fprintf(stderr,"\n");
  632.  
  633.    (*alen) = pos;
  634.  
  635. }
  636.  
  637. static void pdel(short k)
  638. {
  639.         if(last_print<0)
  640.                 last_print = displ[print_ptr-1] -= k;
  641.         else
  642.                 last_print = displ[print_ptr++] = -(k);
  643. }
  644.  
  645. static void padd(short k)
  646. {
  647.  
  648.         if(last_print<0) {
  649.                 displ[print_ptr-1] = k;
  650.                 displ[print_ptr++] = last_print;
  651.         }
  652.         else
  653.                 last_print = displ[print_ptr++] = k;
  654. }
  655.  
  656. static void palign(void)
  657. {
  658.         displ[print_ptr++] = last_print = 0;
  659. }
  660.  
  661.  
  662. /* static int HH[MAXLEN+1], DD[MAXLEN+1], RR[MAXLEN+1], SS[MAXLEN+1]; */
  663. extern int * HH, * DD, * RR, * SS;
  664.  
  665. static int pdiff(short A,short B,short M,short N)
  666. {
  667.         short midi,midj,type;
  668.         int midh;
  669.  
  670.         static int t, tl, g, h;
  671.  
  672. {    static short i,j;
  673.         static int hh, f, e, s;
  674.  
  675. /* Boundary cases: M <= 1 or N == 0 */
  676. if (debug>2) fprintf(stderr,"A %d B %d M %d N %d midi %d\n", 
  677. (pint)A,(pint)B,(pint)M,(pint)N,(pint)M/2);
  678.  
  679. /* if sequence B is empty....                                            */
  680.  
  681.         if(N<=0)  {
  682.  
  683. /* if sequence A is not empty....                                        */
  684.  
  685.                 if(M>0) {
  686.  
  687. /* delete residues A[1] to A[M]                                          */
  688.  
  689.                         pdel(M);
  690.                 }
  691.  
  692.                 return(-gap_penalty1(A,B,M));
  693.         }
  694.  
  695. /* if sequence A is empty....                                            */
  696.  
  697.         if(M<=1) {
  698.                 if(M<=0) {
  699.  
  700. /* insert residues B[1] to B[N]                                          */
  701.  
  702.                         padd(N);
  703.                         return(-gap_penalty2(A,B,N));
  704.                 }
  705.  
  706. /* if sequence A has just one residue....                                */
  707.  
  708.                 midh =  -gap_penalty2(A+1,B,1)-gap_penalty1(A+1,B+1,N);
  709.                 hh   =  -gap_penalty1(A,B+1,N)-gap_penalty2(A+1,B+N,1);
  710.                 if(hh>midh) midh = hh;
  711.                 midj = 0;
  712.                 for(j=1;j<=N;j++) {
  713.                         hh = -gap_penalty1(A,B+1,j-1) + prfscore(A+1,B+j)
  714.                             -gap_penalty1(A+1,B+j+1,N-j);
  715.                         if(hh>midh) {
  716.                                 midh = hh;
  717.                                 midj = j;
  718.                         }
  719.                 }
  720.  
  721.                 if(midj==0) {
  722.                         pdel(1);
  723.                         padd(N);
  724.                 }
  725.                 else {
  726.                         if(midj>1) padd(midj-1);
  727.                         palign();
  728.                         if(midj<N) padd(N-midj);
  729.                 }
  730.                 return midh;
  731.         }
  732.  
  733.  
  734. /* Divide sequence A in half: midi */
  735.  
  736.         midi = M / 2;
  737.  
  738. /* In a forward phase, calculate all HH[j] and HH[j] */
  739.  
  740.         HH[0] = 0.0;
  741.         t = -open_penalty1(A,B+1);
  742.         tl = -ext_penalty1(A,B+1);
  743.         for(j=1;j<=N;j++) {
  744.                 HH[j] = t = t+tl;
  745.                 DD[j] = t-open_penalty2(A+1,B+j);
  746.         }
  747.  
  748.         t = -open_penalty2(A+1,B);
  749.         tl = -ext_penalty2(A+1,B);
  750.         for(i=1;i<=midi;i++) {
  751.                 s = HH[0];
  752.                 HH[0] = hh = t = t+tl;
  753.                 f = t-open_penalty1(A+i,B+1);
  754.  
  755.                 for(j=1;j<=N;j++) {
  756.                     g = open_penalty1(A+i,B+j);
  757.                     h = ext_penalty1(A+i,B+j);
  758.                         if ((hh=hh-g-h) > (f=f-h)) f=hh;
  759.                     g = open_penalty2(A+i,B+j);
  760.                     h = ext_penalty2(A+i,B+j);
  761.                         if ((hh=HH[j]-g-h) > (e=DD[j]-h)) e=hh;
  762.                         hh = s + prfscore(A+i, B+j);
  763.                         if (f>hh) hh = f;
  764.                         if (e>hh) hh = e;
  765.  
  766.                         s = HH[j];
  767.                         HH[j] = hh;
  768.                         DD[j] = e;
  769.  
  770.                 }
  771.         }
  772.  
  773.         DD[0]=HH[0];
  774.  
  775. /* In a reverse phase, calculate all RR[j] and SS[j] */
  776.  
  777.         RR[N]=0.0;
  778.         tl = 0.0;
  779.         for(j=N-1;j>=0;j--) {
  780.                 g = -open_penalty1(A+M,B+j+1);
  781.                 tl -= ext_penalty1(A+M,B+j+1);
  782.                 RR[j] = g+tl;
  783.                 SS[j] = RR[j]-open_penalty2(A+M,B+j);
  784.         }
  785.  
  786.         tl = 0.0;
  787.         for(i=M-1;i>=midi;i--) {
  788.                 s = RR[N];
  789.                 g = -open_penalty2(A+i+1,B+N);
  790.                 tl -= ext_penalty2(A+i+1,B+N);
  791.                 RR[N] = hh = g+tl;
  792.                 t = open_penalty1(A+i,B+N);
  793.                 f = RR[N]-t;
  794.  
  795.                 for(j=N-1;j>=0;j--) {
  796.                     g = open_penalty1(A+i,B+j+1);
  797.                     h = ext_penalty1(A+i,B+j+1);
  798.                         if ((hh=hh-g-h) > (f=f-h-g+t)) f=hh;
  799.                         t = g;
  800.                     g = open_penalty2(A+i+1,B+j);
  801.                     h = ext_penalty2(A+i+1,B+j);
  802.                         hh=RR[j]-g-h;
  803.                         if (i==(M-1)) e=SS[j]-h;
  804.                         else e=SS[j]-h-g+open_penalty2(A+i+2,B+j);
  805.                         if (hh > e) e=hh;
  806.                         hh = s + prfscore(A+i+1, B+j+1);
  807.                         if (f>hh) hh = f;
  808.                         if (e>hh) hh = e;
  809.  
  810.                         s = RR[j];
  811.                         RR[j] = hh;
  812.                         SS[j] = e;
  813.  
  814.                 }
  815.         }
  816.         SS[N]=RR[N];
  817.  
  818. /* find midj, such that HH[j]+RR[j] or DD[j]+SS[j]+gap is the maximum */
  819.  
  820.         midh=HH[0]+RR[0];
  821.         midj=0;
  822.         type=1;
  823.         for(j=0;j<=N;j++) {
  824.                 hh = HH[j] + RR[j];
  825.                 if(hh>=midh)
  826.                         if(hh>midh || HH[j]!=DD[j] && RR[j]==SS[j]) {
  827.                                 midh=hh;
  828.                                 midj=j;
  829.                         }
  830.         }
  831.  
  832.         for(j=N;j>=0;j--) {
  833.                 g = open_penalty2(A+midi+1,B+j);
  834.                 hh = DD[j] + SS[j] + g;
  835.                 if(hh>midh) {
  836.                         midh=hh;
  837.                         midj=j;
  838.                         type=2;
  839.                 }
  840.         }
  841. }
  842.  
  843. /* Conquer recursively around midpoint                                   */
  844.  
  845.  
  846.         if(type==1) {             /* Type 1 gaps  */
  847. if (debug>2) fprintf(stderr,"Type 1,1: midj %d\n",(pint)midj);
  848.                 pdiff(A,B,midi,midj);
  849. if (debug>2) fprintf(stderr,"Type 1,2: midj %d\n",(pint)midj);
  850.                 pdiff(A+midi,B+midj,M-midi,N-midj);
  851.         }
  852.         else {
  853. if (debug>2) fprintf(stderr,"setting to 0: midj %d %d\n",(pint)B+midj,(pint) B+midj+1);
  854.                 profile2[B+midj][GAPCOL] = 0.0;
  855. if (debug>2) fprintf(stderr,"Type 2,1: midj %d\n",(pint)midj);
  856.                 pdiff(A,B,midi-1,midj);
  857.                 pdel(2);
  858. if (debug>2) fprintf(stderr,"Type 2,2: midj %d\n",(pint)midj);
  859.                 pdiff(A+midi+1,B+midj,M-midi-1,N-midj);
  860.         }
  861.  
  862.         return midh;       /* Return the score of the best alignment */
  863. }
  864.  
  865. /* calculate the score for opening a gap at residues A[i] and B[j]       */
  866.  
  867. static int open_penalty1(short i, short j)
  868. {
  869.    int g;
  870.  
  871.    if (i==0 || i==prf_length1) return(0);
  872.  
  873.    g = profile2[j][GAPCOL] + profile1[i][GAPCOL];
  874.    return(g);
  875. }
  876.  
  877. /* calculate the score for extending an existing gap at A[i] and B[j]    */
  878.  
  879. static int ext_penalty1(short i, short j)
  880. {
  881.    int h;
  882.  
  883.    if (i==0 || i==prf_length1) return(0);
  884.  
  885.    h = profile2[j][LENCOL];
  886.    return(h);
  887. }
  888.  
  889. /* calculate the score for a gap of length k, at residues A[i] and B[j]  */
  890.  
  891. static int gap_penalty1(short i, short j, short k)
  892. {
  893.    short ix;
  894.    int gp;
  895.    int g, h;
  896.  
  897.    if (k <= 0) return(0);
  898.    if (i==0 || i==prf_length1) return(0);
  899.  
  900.    g = profile2[j][GAPCOL] + profile1[i][GAPCOL];
  901.    for (ix=0;ix<k && ix+j<prf_length2;ix++)
  902.       h = profile2[ix+j][LENCOL];
  903.  
  904.    gp = g + h * k;
  905.    return(gp);
  906. }
  907. /* calculate the score for opening a gap at residues A[i] and B[j]       */
  908.  
  909. static int open_penalty2(short i, short j)
  910. {
  911.    int g;
  912.  
  913.    if (j==0 || j==prf_length2) return(0);
  914.  
  915.    g = profile1[i][GAPCOL] + profile2[j][GAPCOL];
  916.    return(g);
  917. }
  918.  
  919. /* calculate the score for extending an existing gap at A[i] and B[j]    */
  920.  
  921. static int ext_penalty2(short i, short j)
  922. {
  923.    int h;
  924.  
  925.    if (j==0 || j==prf_length2) return(0);
  926.  
  927.    h = profile1[i][LENCOL];
  928.    return(h);
  929. }
  930.  
  931. /* calculate the score for a gap of length k, at residues A[i] and B[j]  */
  932.  
  933. static int gap_penalty2(short i, short j, short k)
  934. {
  935.    short ix;
  936.    int gp;
  937.    int g, h;
  938.  
  939.    if (k <= 0) return(0);
  940.    if (j==0 || j==prf_length2) return(0);
  941.  
  942.    g = profile1[i][GAPCOL] + profile2[j][GAPCOL];
  943.    for (ix=0;ix<k && ix+i<prf_length1;ix++)
  944.       h = profile1[ix+i][LENCOL];
  945.  
  946.    gp = g + h * k;
  947.    return(gp);
  948. }
  949.